OCT  1  1981 

OCT  1  3  1981 


AEDC-TR  80  56 

c-'b 


wm 

IP 


Computation  of  Three-Dimensional 
Transonic  Internal  Flow 
in  Cylindrical  Coordinates 


W.  J.  Phares 
ARO,  Inc. 


August  1981 


Final  Report  for  Period  September  30,  1979  -  October  1,  1980 


Approved  for  public  release;  distribution  unlimited. 


Pi'OpSi'V 

dBHARY 

F406G0-6t-C-OtlO4 


ARNOLD  ENGINEERING  DEVELOPMENT  CENTER 
ARNOLD  AIR  FORCE  STATION,  TENNESSEE 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 


NOTICES 


When  U.  S.  Government  drawings,  specifications,  or  other  data  are  used  for  any  purpose  other 
than  a  definitely  related  Government  procurement  operation,  the  Government  thereby  incurs  no 
responsibility  nor  any  obligation  whatsoever,  and  the  fact  that  the  Government  may  have 
formulated,  furnished,  or  in  any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is 
not  to  be  regarded  by  implication  or  otherwise,  or  in  any  manner  licensing  the  holder  or  any 
other  person  or  corporation,  or  conveying  any  rights  or  permission  to  manufacture,  use,  or  sell 
any  patented  invention  that  may  in  any  way  be  related  thereto. 


Qualified  users  may  obtain  copies  of  this  report  from  the  Defense  Technical  Information  Center. 


References  to  named  commercial  products  in  this  report  are  not  to  be  considered  in  any  sense 
as  an  indorsement  of  the  product  by  the  United  States  Air  Force  or  the  Government. 


This  report  has  been  reviewed  by  the  Office  of  Public  Affairs  (PA)  and  is  releasable  to  the 
National  Technical  Information  Service  (NTIS).  At  NTIS,  it  will  be  available  to  the  general 
public,  including  foreign  nations. 


APPROVAL  STATEMENT 

This  report  has  been  reviewed  and  approved. 


KEITH  L.  KUSHMAN 


Directorate  of  Technolpgy 
Deputy  for  Operations 


Approved  for  publication: 
FOR  THE  COMMANDER 


MARION  L.  LASTER 
Director  of  Technology 
Deputy  for  Operations 


REPORT  DOCUMENTATION  PAGE 


HEAD  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


2  GOVT  ACCESSION  NO.  3  RECIPIENT’S  CATALOG  NUMBER 


1  REPORT  NUMBER 

AEDC-TR-80-56 


*  T'TLE  ( tnd  Sutlillw)  S  TVPE  OF  REPORT  A  PERIOD  COVERED 

COMPUTATION  OF  THREE-DIMENSIONAL  TRANSONIC  Final  -  Sept.  30,  1979 


INTERNAL  FLOW  IN  CYLINDRICAL  COORDINATES 


f  author^ 

W.  J.  Phares,  ARO,  Inc.,  a  Sverdrup 
Corporation  Company 


-  Oct.  1,  1980 


6  PERFORMING  ORG.  REPORT  NUMBER 


B  CONTRACT  OR  GRANT  NUMBERfsJ 


10  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Program  Element  65807F 


<2.  REPORT  DATE 


August  1981 


9  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS  ' 

Arnold  Engineering  Development  Center/DOT 

Air  Force  Systems  Command 

Arnold  Air  Force  Station,  Tennessee  37389 


CONTROLLING  OFFICE  NAME  AND  ADORESS 

Arnold  Engineering  Development  Center/DOS 
Air  Force  Systems  Command 
Arnold  Air  Force  Station,  Tennessee  37389  33 

>«  MONITORING  AGENCY  NAME  »  AODRESSf/f  dlltmrvnt  from  Contrail  In,  Ot(lct)  15  SECURITY  CLA5S.  <al  Ihl,  rspo/rj 

UNCLASSIFIED 

ISA.  DECLASSIFICATION -’DOWN  GRADING 
SCHEDULE  ..  .  . 

_ _ _ _ _  N/A 

T*  OlSTPlBuTlON  STATEMENT  fol  Ihlr  Rtpail)  " - - - - 

Approved  for  public  release;  distribution  unlimited. 


1?  DISTRIBUTION  statement  (ol  the  mbmtrmct  tmfermd  in  Block  20,  if  different  from  Rtpait) 


is  supplementary  notes 
Available  in  Defense  Technical  Information  Center  (DTIC) . 


■5  KEY  WORDS  (Continue  on  r avetee  etde  tt  nmcm 

computer  programs 
transonic  flow 
exhaust  nozzles 
three-dimensional  flow 


aery  and  Identify  by  block  number) 

propulsion  systems 
test  methods 


20  ABSTRACT  (Contlrmm  or  imvotte  tldt  II  ntcetl mrr  and  Idanllly  by  Hoc*  number)  ' 

The  Euler  equations  for  the  inviscid  rotational  flow  of  a 
perfect  gas  are  solved  for  three-dimensional  subsonic-transonic 
internal  flow  by  use  of  a  time-dependent  numerical  technique.  The 
numerical  approach,  which  is  an  extension  and  modification  of  that 
of  Cline  for  two-dimensional  flows,  is  basea  on  the  use  of  the 
MacCormack  finite-difference  method  for  the  interior  field 
points.  The  reference-plane  method  of  characteristics  is  used  for 


DO  i  jan*73  1473 


EDITION  OF  I  NOV  55  IS  OBSOLETE 


UNCLASSIFIED 


UNCLASSIFIED 


20.  ABSTRACT  (Continued) 

coupling  the  interior  field  solution  to  the  boundary  points. 
Both  the  basic  equations  and  the  numerical  procedures  are 
described,  as  is  the  computer  program  which  was  written  in 
FORTRAN  IV  language  for  either  the  Cray-1  or  the  IBM  370/165 
computer.  As  presently  written,  the  program  is  applicable 
to  the  computation  of  three-dimensional  flow  in  both  axi- 
symmetric  and  relatively  simple  three-dimensional  nozzle 
geometries.  The  validity  of  the  computer  program  was  estab¬ 
lished  by  computing,  in  various  ways,  an  axisymmetric  nozzle 
flow  as  a  three-dimensional  flow;  the  numerical  results  are 
in  good  agreement  with  the  results  from  a  well-established 
computer  program  for  axisymmetric  flow. 


APB  C 

Ai-cJC  APS  T*nn 


UNCLASSIFIED 


A  ED  C-TR  -80-56 


PREFACE 

The  work  reported  herein  was  conducted  by  the  Arnold  Engineering  Development 
Center  (AEDC),  Air  Force  Systems  Command  (AFSC).  The  results  presented  were  obtained 
by  ARO,  Inc.,  AEDC  Group  (a  Sverdrup  Corporation  Company),  operating  contractor  for 
the  AEDC,  AFSC,  Arnold  Air  Force  Station,  Tennessee.  Elton  R.  Thompson  was  the  Air 
Force  project  manager.  The  work  was  done  under  ARO  Project  No.  E32A-01,  and  the 
manuscript  was  submitted  for  publication  October  6,  1980. 
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1.0  INTRODUCTION 

Theoretical  calculations  of  exhaust  nozzle  performance  are  often  required  to  aid  in  the 
evaluation  of  propulsion  systems  tested  at  the  Arnold  Engineering  Development  Center 
(AEDC).  Computer  codes  have  been  available  for  axisymmetric  transonic  nozzle  flow  (Refs. 
1  and  2}  and  the  supersonic  portion  of  axisymmetric  flow  fields  can  be  calculated  with  the 
well-known  method  of  characteristics. 

Even  with  axisymmetric  nozzle  geometries,  the  exhaust  nozzle  flow  field  is  often  three- 
dimensional  (3-D)  because  of  asymmetries  in  the  entrance  flow;  future  tests  at  AEDC  will 
involve  engines  with  truly  three-dimensional  nozzle  geometry.  Consequently,  there  is  a 
pressing  need  for  the  capability  to  make  computations  of  three-dimensional  nozzle  flow 
fields.  A  computer  program  has  been  developed  for  calculating  the  supersonic  flow  in  3-D 
nozzles  (Ref.  3),  but  no  program  has  been  available  for  computing  the  subsonic-transonic 
flow  in  such  nozzles. 

In  the  present  study,  a  computer  program  has  been  developed  for  the  inviscid  subsonic- 
transonic  flow  in  three-dimensional  propulsion  nozzles.  Of  course,  it  is  not  limited  to 
computation  of  propulsion  nozzle  performance;  it  is  also  applicable  to  many  other  3-D 
internal  flows.  Because  of  the  unavailability  of  detailed  experimental  data  or  other  analytical 
solutions  for  3-D  nozzles,  the  current  program  was  evaluated  using  axisymmetric  nozzle 
geometry  with  the  computational  axis  offset  from  the  nozzle  axis  to  simulate  3-D  flow  fields. 


2.0  ANALYSIS 


2.1  APPROACH 

The  Eulerian  equations  in  nonconservative  form  are  solved  for  the  three-dimensional, 
inviscid  rotational  flow  of  a  perfect  gas  with  a  time-dependent  numerical  technique.  The 
objective  of  this  research  was  to  extend  and  modify  the  successful  two-dimensional  (2-D) 
Cline  method  (Ref.  1)  to  provide  a  3-D  capability. 

According  to  Cline,  long  computation  times  associated  with  time-dependent  2-D 
calculations  are  usually  required  because  inefficient  algorithms  or  poor  treatment  of 
boundaries  demand  excessively  fine  computational  meshes.  Cline,  using  the  MacCormack 
scheme  coupled  with  characteristic  boundary  conditions,  produced  a  2-D  code  with 
reasonable  computational  times.  His  success  motivated  the  present  approach  to  the  3-D 
problem.  An  outline  of  the  approach  follows: 
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1  *  The  Eulerian  equations  in  nonconservative  form  are  solved. 

2.  Interior  mesh  point  properties  are  computed  using  the  efficient  MacCormack 
finite  difference  scheme. 

3.  The  inlet  and  wall  boundary  mesh  point  properties  are  calculated  using  a 
reference-plane  characteristic  technique. 

4.  Exit  mesh  point  properties  are  calculated  using  linear  extrapolation  for 
supersonic  flow  and  a  characteristics  scheme  for  the  subsonic  case. 

5.  Singularities  along  the  z-axis,  attributable  to  the  choice  of  a  cylindrical 
coordinate  system,  are  avoided  by  excluding  the  z-axis  from  the  flow  field 
calculations. 

6.  Accumulation  of  truncation  error  is  reduced  by  an  alternating  scheme  for  the 
backward-forward  option  in  MacCormack’s  method. 

Note  that  items  4  and  5  are  different  from  Cline’s  approach.  The  physical  (z,  r,  0)  space  is 
transformed  into  a  right-circular  cylinder  computational  domain  (£,  t),  f)  by  a  coordinate 
transformation.  The  computational  mesh  is  uniform  in  the  £,  r\,  and  f  directions. 

The  physical  space  mesh  has  equal  spacing'  in  the  axial  (z)  and  circumferential  (0) 
directions,  whereas  that  in  the  radial  (r)  directions  may  be  unequal. 

The  temporal  step  size,  At,  is  controlled  by  the  Courant-Fredricks-Lewy  (CFL) 

condition.  In  general,  no  smoothing  or  damping  techniques  are  required  to  maintain 
stability. 


2.2  GOVERNING  EQUATIONS 


The  governing  equations  in  cylindrical  coordinates  for  time-dependent  inviscid, 
adiabatic  3-D  flow  of  an  ideal  gas  are 


Continuity: 


Momentum: 


Di 


+ 


V  »  (p  v  )  =  0 


_1  <?p 

P  dr 


0 


(1) 


(2) 
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[D>*  ♦'■'•l,1  n  ,o 

L  ot  r  J  Pr  de 

(3) 

l  +  =  ° 

L  Di  J  p  dz 

(4) 

Energy: 

W-.*  ^  =  o 

LDiJ  Dt 

(5) 

State: 

p  =  pRT 

(6) 

where 

DO  dO  do  '0  d()  dO 

=  -  1  V  -  ^ - +  v  - 

Dt  di  r  dr  r  S6  x  dz 

(7) 

and 

v  •  <pv)  =  r  “  (v  vr)  + 1 ~  (pvff)  + 1-  (p  vB) 


(8) 


and  where  p  is  the  density,  is  the  axial  velocity,  vr  is  the  radial  velocity,  vfl  is  the 
circumferential  velocity,  p  is  the  pressure,  a  is  the  local  speed  of  sound,  t  is  the  time,  and  z,  r, 
0  are  the  axial,  radial,  and  circumferential  coordinates. 

The  physical  space  (z,  r,  0)  is  mapped  into  a  right-circular  cylinder  computation  domain 
(£,  tj,  vf)  by  the  following  coordinate  transformation: 


f  -  * 


where  rw  denotes  the  wall  values.  In  the  (£,  v,  f)  coordinate  system,  Eqs.  (1)  through  (5) 
become  (where  the  subscripts  outside  parentheses  indicate  partial  differentiation): 
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Continuity: 


Momentum: 


Energy: 


8vs 


+  V(p)/?  +  r  V0  ( p )Jf  +  —  ( p +  vz(p)^ 


°2P 


+  p(vz)^  -«lP(v^TP^  =  0 


(9) 


(Vr),  *V(V-),  +  Vz(Vr)f*ST  (V>)f  *  f  (PV 


*? 


CL  r)  V/n 

w  *  r?  /  \  n 

- +  - (  V  ^  =  o 

r  r  V  >V 


(10) 


(vfl)t’v(vfl),  +  -7^-  (V0)/“T"  (ve)^+vzCv0)^ 


\  '$  fl2  |  v  5  ,  .  . 

+• -  +  —  (p)„  +  —  tpV  =  o 

r  (Or  r  7  pr  r  <, 


(ID 


ve 


5'v 


(’.),* v  (••  V  -  Vs-  ('• ),  *  -r-(’> V 0.)f 


%  <pif4  w,  - 0 


(12) 


{p)j  t-  V  (p)^  +  V0  (r)^  +  7  V0  (p)^+  Vz  (p)£ 


+  v 


[(p\  + 

(p)fj=  0 


—  a2vf)  &  vfl 

v  +  (p),  +  — r^-  (p)c 


(13) 
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where 

V  =  /3  vr  +  Oj  Vz 

0  -  (i?)r  -  p- 

W 

“i  =  Wz  =  ~fa(r»\ 

a2  =  We  =  -^(rw)0 

and 

(fl,  -  (Or  "0,8-  (i)S  • 

3.0  NUMERICAL  METHOD 

The  computational  space  is  divided  into  four  types  of  mesh  points:  interior,  inlet,  exit, 
and  wall. 

3.1  INTERIOR  MESH  POINTS 

The  interior  mesh  point  values  are  calculated  by  the  second-order  accurate,  explicit 
MacCormack  finite-difference  method,  using  an  alternating  backward-forward  scheme.  The 
alternating  backward- forward  scheme  used  is  illustrated  in  Fig.  1.  The  governing  equations 
are  left  in  nonconservative  form. 

3.2  INLET  MESH  POINTS 

The  inlet  mesh  points  for  subsonic  flow  are  computed  using  a  second-order,  reference 
plane  characteristic  scheme.  In  this  method,  the  partial  derivatives  with  respect  to  rj,  f  are 
computed  in  the  initial  value  and  solution  surfaces  using  noncentered  differences  as  in  the 
MacCormack  scheme.  These  approximations  are  then  transferred  to  the  right-hand  side 
(RHS)  in  the  governing  equations,  and  the  resulting  system  of  equations  is  solved  in  the 
constant  tj,  f  reference  planes  using  a  two-in  dependent-variable  characteristics  scheme. 

The  characteristics  relations  derived  in  Appendix  A  that  relate  the  interior  flow  to  the 
nozzle  inlet  flow  are 

dp-pa<lvz  =  (a2  ifjl -pa 

for  !r  =  l'z  “  a  (14) 
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where  ^  is  defined  in  Appendix  A. 
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1  -  First  Step  of  Two-Step  Method 

2  -  Second  Step  of  Two-Step  Method 
B  -'Backward  Difference 

F  -  Forward  Difference 

Figure  1.  Alternating  scheme  for  the  backward-forward 
option  in  MacCormack's  method. 

Use  of  a  reference-plane  characteristics  scheme  requires  the  specification  of  the  inlet  flow 
angles  as  well  as  the  stagnation  conditions.  The  equations  relating  the  total  and  static 
conditions  are 


pt/p  . 


Tr/T  =■  ]  +(y-i)  M2/  2  (16) 

where  y  is  the  ratio  of  specific  heats,  M  is  the  Mach  number,  T  is  the  temperature,  and  the 
subscript  T  denotes  the  total  (stagnation)  condition. 

Equation  (14)  is  solved  using  standard  characteristic  techniques,  with  the  t  terms 
evaluated  in  the  initial  value  plane.  Equations  (14),  (15),  and  (16),  along  with  the  inlet  flow 

angles  and  the  equation  of  state,  form  a  system  of  five  equations  for  the  five  dependent 
variables. 
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3.3  EXIT  MESH  POINTS 


For  subsonic  exit  flow,  a  reference-plane  characteristic  scheme  similar  to  the  inlet  scheme 
is  used.  The  exit  pressure  is  specified.  The  characteristic  relations  relating  the  interior  flow  to 
the  nozzle  exit  flow,  derived  in  Appendix  A,  are 


and 


d  vr  =  iA2  dt 

d  *0  =  1A3  dl 

dp  —  a2  dp  =  dt 


for 


dt 


(17) 

(18) 
(19) 


(20) 


Equations  (17),  (18),  (19),  and  (20),  along  with  the  exit  pressure  condition,  form  a  system  of 
five  equations  for  the  dependent  variables. 


For  supersonic  flow,  the  flow  conditions  at  the  exit  mesh  points  are  computed  by  linear 
extrapolation. 

s 

3.4  WALL  MESH  POINTS 


The  wall  mesh  points  are  also  computed  using  a  reference-plane  characteristic  scheme.  In 
this  scheme,  the  derivatives  with  respect  to  £,  f  are  approximated  and  the  resulting  system  of 
equations  is  solved  in  the  £,  f  =  constant  reference  planes. 


The  characteristic  relations  derived  in  Appendix  B  which  relate  the  interior  flow  to  the 
flow  at  the  nozzle  walls  are 


dp  —  a2  dp  = 

dvr  -  ~  dvfl  = 
dv.  — —  dv_  = 


*^5  d*- 


**)dl 


dl  (  for  £7  _  y 

dt 


(21) 

(22) 

(23) 
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.  a]  ,  ■  B  a, 

dp  +  pa  —  dv  +  pa  -C.  dv  +  pa  —  dvfl 
a*  a*  r  a*  “ 

=  i  a2  +pa  -£  +  «fr3  +  pa  fi 


where 


and 


for  *1  =  V  + 

dt 


(24) 


(p+r'-ij 


V  =  J8vr 


«2  v0 


+  “L  v« 


Equations  (21),  (22),  (23),  and  (24),  along  with  the  wall  boundary  condition,  form  a  system 
of  five  equations  for  the  dependent  variables.  The  wall  boundary  condition  is  given  by 

V,F.  +  Vrl’r^^  =  0 

Where  F  =  constant  defines  the  wall  suface  and  Fz,  Fr,  and  F,  are  the  corresponding  partial 
derivatives. 

3,5  TIME  STEP  SIZE 

The  time  step  size.  At,  is  controlled  by  the  well-known  CFL  condition  which  can  be 
expressed  as 


At  <  \j  {(  V  +  a)  [  ]  /(Atj/jS) 2  -  1  /( A0 2  +  1  A\0 2]  }  (25; 

where  V  is  the  velocity  magnitude.  Experience  gained  in  the  present  study  indicates  that 
values  of  A  from  1.0  to  2.0  are  satisfactory. 


4.0  RESULTS  AND  DISCUSSION 

Several  adiabatic  flow  cases  were  selected  to  verify  the  numerical  approach  described  in 
the  preceding  sections.  The  test  cases  were: 
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1.  Steady-state  uniform  flow  in  a  cylindrical  duct  with  an  initial  disturbance  at  one 
mesh  point. 

2.  Axisymmetric  steady  flow  in  a  choked  converging-diverging  nozzle.  The  nozzle 
geometry  was  the  same  as  for  Cline’s  Case  1  (Ref.  2)  and  consisted  of  a  45-deg 
conical  inlet  section,  a  circular  arc  throat  section,  and  a  15-deg  conical  diverging 
section  (Fig.  2). 

3.  Case  2  with  a  nonaxisymmetric  swirl  induced  by  a  “bump”  in  the  otherwise 
uniform  inlet  total  pressure  distribution. 

4.  Case  2  with  the  computational  axis  offset  various  amounts  from  the  nozzle  axis, 
which  simulates  complex  3-D  flow  fields. 


Figure  2.  Axisymmetric  nozzle  geometry. 

Cases  I  and  2  were  used  to  debug  and  verify  the  computer  program  for  rather  simple 
flows  with  known  steady-state  solutions.  For  Case  1,  the  solution  relaxed  smoothly  toward 
uniform  flow.  (In  none  of  the  test  cases  did  the  solution  at  very  large  times  tend  to  diverge 
from  the  known  steady  solution).  The  steady-state  solution  for  Case  2  was  negligibly 
different  from  the  results  reported  by  Cline. 


13 


AEDC-TR-80-56 


A  converged  steady  solution  obtained  for  Case  3  indicates  that  the  numerical  approach  is 
applicable  to  predicting  the  effeci  of  nonsymmetrical  inlet  profiles  on  the  flow  in 
axisymmetric  nozzles. 


Figure  3.  Computational  mesh  in  throat  plane  for 
case  4  {20-percent  offset). 


For  Case  4,  the  computational  axis  was  offset  from  the  axisymmetric  nozzle  axis  by  6, 
10,  and  20  percent  of  the  nozzle  throat  radius,  r*.  The  computational  mesh  in  the  throat 
plane  is  illustrated  in  Fig.  3  for  20-percent  offset.  The  meshes  used  and  the  run  times  on  the 
Cray-1  computer  are  summarized  in  Fig.  4  for  the  three  axis  offsets.  For  all  three  offsets,  the 
steady-state  static  pressures  throughout  the  flow  field  differed  less  than  one  percent  from 
those  calculated  with  the  axisymmetric  Cline  program  (which  correlate  well  with 
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experiment).  These  computations  of  axisymmetric  nozzle  flow,  done  the  “hard  way”  with 
the  3-D  program,  indicate  that  the  program  is  acceptably  accurate  and  is  capable  of 
computing  3-D  subsonic-transonic  flow. 

With  one  exception,  all  solutions  were  obtained  without  application  of  smoothing  or 
damping  techniques  to  maintain  stability.  The  exception  was  Case  4c,  with  20-percent  axis 
offset;  in  that  one  case,  stability  could  not  be  achieved  without  a  modification  of 
MacCormack’s  scheme  which  adds  a  second  order  truncation  error  and  has  the  effect  of  an 
artificial  viscosity.  Whether  the  use  of  such  damping  techniques  can  be  avoided  by  proper 
mesh  selection  remains  a  subject  for  further  investigation. 


Run  Time, 


Case 

Of fset /r* 

Grid 

min  (Cray) 

4a 

0.06 

21  x  8  x 

8 

0.  37 

4b 

0.10 

21  x  12  x 

12 

3. 70 

4c 

0.20 

21  x  12  x 

18 

10.1 

Figure  4.  Axisymmetric  nozzle  with  offset  computational  axis. 


5.0  CONCLUDING  REMARKS 

The  Cline  numerical  approach,  which  is  applicable  to  planar  and  axisymmetric  nozzle 
flows,  has  been  successfully  extended  to  three-dimensional  flows  in  this  study.  Results 
obtained  for  the  various  test  flows  indicate  that  the  present  computer  program  is  applicable 
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to  prediction  of  the  effect  of  nonsymmetric  inlet  profiles  on  the  flow  in  axisymmetric  nozzle 
geometries.  In  addition,  the  computer  program  seems  well-suited  to  the  computation  of  the 
flow  fields  in  relatively  simple  3-D  nozzle  geometries.  However,  additional  work  remains  to 
be  done  in  extending  the  program  to  arbitrary  and  complex  3-D  geometries. 
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APPENDIX  A 

CHARACTERISTIC  RELATIONS  FOR  INLET  AND  EXIT  MESH  POINTS 

Governing  Equations 

The  governing  Eqs.  (18)  through  (22)  can  be  written  as 

H  +  V  (p )r}  +  &p  (v^  +  (v0)  +  a  ,  p  (vz  )  = 


"  {~r  {p](  +  + ~r  p  (ve)y 4 p  (v0  * 4  V 


(A-l) 


(A-2) 


(v0)t4  v(ve)I}  +  ~  =  -{~(ve)c+vz(v0)f+-7^  + 


Pr 


(p)^ 


(A-3) 


WtlV(v0,7  {p)7 =  'I 


{  8  ve(\)c 


T -  +  vz(vz  +  p  (P)f 


(A-4) 


(p)i+V(p)7-a2((P)l+V(^)  =  J 


-  \Q  (p>£  +  Vj,  <p>£ 


S  V/ 


—  (P^+Vz(p)^r 


(A-5) 


where  V  —  0vr  +  a2  vj/r  +  a  i ,  vz  is  the  tangential  condition 


Defining 


=  RHS  of  (A-l) 
2  =  RHS  of  (A-2) 
=  RHS  of  (A-3) 


=  RHS  of  (A-4) 


=  RHS  of  (A-5) 
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then  the  system  of  governing  Equations  {A-l)  through  (A-5)  can  be  written  in  vector 
notation  as 


where 


AIjtB  W'jj  -  V 


p 

^1 

vr 

-* 

** 

vd 

and  F  = 

^3 

v* 

^4 

P 

^5 

(A-6) 


and  with  the  A  and  B  coefficient  matricies  as 


A 


1  0  0  0  0 

0  1  0  '  0  0 

0  0  10  0 

0  0  0  10 


-a2  0001 


and 


B 


V 

0 

0 


ft> 


“2  P 

r 

v  o 


0  V 


-a3  V 


0  0 


a\P  0 
0  ft/p 


Pr 

P 

V 
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Defining  Characteristics 


The  theory  of  eigenvalues  is  used  to  derive  the  equations  of  the  characteristics.  The 
eigenvalues  are  defined  as 


X  =  *L 

dt 


(A- 7) 


The  total  derivative  of  the  W  vector  is 


—  =  w.  -  w„  =  if 


dt 


dt  7 


,  nr  A  W„ 
1  7) 


which,  coupled  with  the  vector  equation,  Eq.  (A-6),  yields 


(B-AA)S  =  F  -  A  — 
7  dt 


(A-8) 


(A-9) 


The  eigenvalues  and  eigenvectors  (A  j)  satisfy  the  inner  product  relation,  based  on  Eq.  (A-9), 
[<B-AjA>»,,A,]  .  A,  ]  .  0  (A-10) 


or 


[  VB-M)*  A,]  -  [f-Af .  A,]  .  0 
where  (B  -  XA)*  signifies  the  transpose  of  the  (B  -  XA)  matrix. 


(All) 


These  matrices  are 

(B-AA)  = 


(V  -  A)  jg p 

0  (V-A) 

0  0 


i2p 


«1  P 
0 


P/p 


-az(V-A)  0 


(V-A)  0  — 

PT 

P 

(V-A) 


0  (V-A) 

0  0 
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and 


(B-AAf  = 


(V-A) 

Pp 

a2P 

r 

“l  P 
0 


0 

V-A 

0 

0 

P 

P 


0 

0 

V-A 

0 

Pr 


0 

0 

0 

V-A 

°1 


-a2(V  -  A) 
0 

0 

0 

V-A 


The  eigenvalues  are  obtained  from  the  determinate  relationship 


|  (B  -  A  A)  *  |  =  (V-A)5-  n2(V-A)  [V(V  -  A)2  +  (V -A)  j  (V-A)  +  (V-A)afJJ 


The  resulting  eigenvalue  equations,  which  describe  the  characteristics,  are 

(V-A)3  =  0 


and 


(V  -  A2)  -  a2 


j  ta>) 


0 

(A- 12) 

(A-13) 

(A- 14) 


Equation  (A-13)  is  easily  solved  for  Xj,  X2  and  X3  characteristics, 

*1  =  a2  "  A3  =  V  (A- 15) 

Equation  (A-14)  is  solved  for  X4  and  X5  characteristics, 

A4  -  A5  =  V  ±aa*  {A- 16) 

where 
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Compatibility  Relations 

To  derive  the  compatibility  relations  along  the  characteristics,  the  eigenvectors  which 
define  the  characteristics  must  satisfy 

[(B-AjA),  Aj  ]  =  0  (A- 17) 

for  \i,2l3  =  V,  or 


The  resulting  eigenvectors  are 


also,  for  X4  5  =  V  ±  aa* 
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the  eigenvectors  are 


A 


4, 


3 


-  (+  aa*)4 


-  /9  p  (±aa*)3 

-Hj  p  (±aa*)3 

-  (iaa*)2  a*  2 


(A-19) 


The  compatibility  relations  along  the  characteristics  are,  using  the  eigenvectors  and  Ea 
(A-ll), 


for 


F  -A  —  , 


1  dL 
d  v 

*2~17 


dv 


8 


^3 - - 

3  <il 

dv 

,  7. 

V*~ 

i/r,+  a2^-^ 
5  dl  dt 


A, 


0 
0 
0 
0 

1  J 


9 


v.  +  a2  -  111  = 


di  di 
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or 


dp 


/<  [  On  tlva  / 

a*  ±  p/3advr  /a*  t  p  a  -  ±  a2  i(/] 


l  a2pa  /  \ 

i  a  &2.ja*  -  - -  “  ±  O]  pa  ^4/a+ jdt 

Compatibility  relations  for  (£,f)  constant  reference  plane  are 


dp  —  a  ~  dp  -=  \!j-  dr 

)dt 


dt 


and 


d  v„ 


(lv  , 


dp  i  pa  j  a  t  p/3  a -  ±  pa2  l» 


-( 


/  «•> 


,  for  =  V 
rlt 


^■5.  ±  n-  ib  j  ±  p  j3a  'iJ2<a>  -  ~~  P  a  ^3  /a"  ±  a  |  pa  tjj 


for  ^ 
di 


=  V+aa1 


dl 
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APPENDIX  B 

CHARACTERISTIC  RELATIONS  FOR  WALL  MESH  POINTS 


Governing  Equations 


The  governing  Eqs.  (18)  through  (32)  can  be  written  as 


(p) 


t+p(v.)f+v«  W*  =  -{v(p)7,  +  -t^  (pV~T  {p)C 

+  Mvr),*“T  (vfl)/TP  Mc+alP  (B_1) 

CVr  \  +  vz  (vr  )f  =  -  {V  (  Vr  ^  +  ~r  (Vr  "  p  (P)IJ  "  ~T  +  (B*2) 


f-  a2  v(9  vr  vf)  1 

(ve)t+  vz(v0)f  =  -|v(ve)/  (we\  +  —(v9)c+  -T“J 


(B-3) 


a2  S 

-  -  (P)„  -  —  (P)/ 

pr  V  pr  h 


(Vz),+  V,  (%)^+  ^  (P)f  =  -{V(v*),+  ^(Vz)>;- 

+  —  j 

p  >7  J 


71  W<)} 


Sva(vz)/ 


(B-4) 


Vp><r,2(5l'V^7i('V 


(B-5) 


where 


V  =  /3  vr  +  a ,  V . 
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and 


<£,  =  R1IS  of  (B-n 
02  =  RMS  of  (H-2) 

6:i  --  RMS  of  (B-3) 
ii4  -  RIlS  of  CQ--1) 

-  RMS  of  (B-5) 

The  theory  of  eigenvalues,  defined  in  Appendix  A,  is  used  to  derive  the  equations  of  the 
characteristics  and  compatibility  relations. 

Equations  (B-l)  through  (B-5)  in  vector  notation  are 


where 


and 


A  «,  f  B  Wg  =  F 


e 

_ 

1“ 

—I 

P 

V'2 

vr 

ft  -= 

'ft 

an 

A  F 

- 

<//2 

P 

^4 

_l 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

A  = 

0 

0 

1 

0 

0 

0 

0 

0 

] 

0 

-a2 

0 

0 

0 

1 
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Vz  0  0  p  0 

0  Vz  0  0  0 

o  0  vz  0  0 

0  0  0  V  - 

z  P 

-a2vz  0  0  0  vz 


then 


(B-AA) 


( v  7  —  A )  0  0  p  0 

0  ( vz -A)  0  0  0 

0  0  (vz-A)  0  0 


0  0  0  (v  —  A )  - 

P 

-a2(vz-A)  0  0  0  (vz-A) 


(B  -  AA)* 


(vz-A) 


0  0 

o  (vz-A)  0 

0  0  (vz-A) 

poo 

ooo 


0 

0 

0 


(vz  -  A) 


J. 

P 


-  a"(vz  -  A) 


0. 

0 

0 


K- A) 


|(H-AXV|  =  (vz-AP-.i2(vz-A)3  ^  (vz-A)3((vz-A)2-a2  )  = 


Aj  = 


a2  = 


= 


vz  ± 
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The  eigenvectors  and  eigenvalues  satisfy 

[tB  -  AA)*  „  A.]  =  0 


for  X|  -  X2  =  X3  =  vz 

0 

0 

0 

p 

0 

then 


0  0  0  0 

r- 

A1 

0  0  0  0 

A2 

0  0  0  0 

A3 

0  0  0  0 

A4 

0  0  1/p  0 

_AS_ 

-  - 

0 

0 

0 

1 

0 

0 

Al.  2,  3  = 

0 

1 

0 

0 

0 

0 

0 

0 

1 

for  X4  5  =  vr  ±  a 


results 


— 

— 

_ 

-  a 

b 

0 

0 

±a3 

0 

+  a 
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0 

0 

0 
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+  a 
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+  a  Aj  +  a3  A5  =  0 
p  A]  +  a  A  4  =  0 


A1  =  ±a2  A5 


let  As  =  1 

Then  A4  =  ±  ga,  A|  -  a2 
and 


±  Pa 
1 


The  compatibility  relations  along  the  characteristics  are,  using  the  eigenvectors 


di 


_ 

Hi 

dv. 


V>3  “ 


Hi 


di 


,  dvz 

**-ir 


i/f.  +  a 


2 


di 


dP 

<ll 
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(*1- £)*<»(** 

Compatibility  relations  for  (ij,£)  constant  reference  planes  are 


^  \  ,  2  dp  dP 

—  - =0 

dt  1  5  dt  dt 
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'l'  2 


dv^  =  dt 
dP  —  a“  dp  =  ^r-dl 


for,  ^  =  v 
dt  z 


dP  ±  pa  dvz  =  ^a2  ^  ±  pa  +-  i dt  |  for, 


M  _  v  +fl 

dt  z 
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A 
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A  ,B 


Fz,  Fr,  Ffl 
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Pr 
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r 

rw 
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vr 


vz 

V" 

Wt,  W,#  Wf 
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NOMENCLATURE 

Premultiplier  for  the  time  step 

Speed  of  sound 

Coefficient  matrices 

Local  wall  slopes 

Static  pressure 

Total  pressure 

Gas  constant 

Radial  coordinate 

Wall  radius 

Nozzle  throat  radius 

Static  temperature 

Time 

Total  temperature 

Velocity  magnitude 

Circumferential  velocity  component 

Radial  velocity  component 

Axial  velocity  component 

Tangential  velocity  condition 

Partial  derivative  dependent  variable  vectors 

Dependent  variable  vector 

Axial  coordinate 

Time  step  size 
Ratio  of  specific  heats 
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0  Circumferential  coordinate 

A  j  Eigenvectors 

\  Eigenvalues 

£,  ij,  f  Transformed  spatial  coordinates  corresponding  to  (z,  r,  0) 

e  Density 

$  Null  vector 

d/dz,  d/dr,  d/d6  Partial  derivative  operators 

*  Transpose 

•  Dot  product  vector  operator 
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